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We present a generalized equations-of-motion method that efficiently calculates energy spectra and 
matrix elements for algebraic models. The method is applied to a 5-dimensional quartic oscillator 
that exhibits a quantum phase transition between vibrational and rotational phases. For certain 
parameters, lOx 10 matrices give better results than obtained by diagonalising 1000 by 1000 matrices. 
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Equations-of-motion methods P, [3, S Q offer an al- 
ternative to diagonalization of a Hamiltonian to deter- 
mine the properties of a quantal system. We consider 
systems with a Hamiltonian expressed in terms of a Lie 
algebra of observables, g. For a finite-dimensional irre- 
ducible representation (finite irrcp) of g, both approaches 
give the same results. However, they become inequiva- 
Icnt when approximations are necessary. For example, 
diagonalizing the Hamiltonian for a system of n coupled 
harmonic oscillators, in a basis of uncoupled harmonic 
oscillator states, gives results that become precise (to 
within computational errors) in an infinite limit whereas 
an equations-of-motion approach, as given by the stan- 
dard random phase approximation (RPA) , is already pre- 
cise in a 2n-dimcnsional space. 

Problems of interest are ones for which an expan- 
sion of the low-energy eigenstates of the Hamiltonian, 
in a given ordered basis for the Hilbert space, converges 
slowly. For example, the low-energy states of strongly- 
deformed rotational nuclei are dominated by components 
from higher shells in contrast to the low-energy states 
of a near-spherical nucleus which may be dominated by 
the valence-shell states of a spherical harmonic-oscillator 
shell-model basis Q. In such a situation, diagonaliza- 
tion of the Hamiltonian in a spherical shell-model basis 
is unlikely to give reliable results. 

The proposed approach avoids the preliminary stage 
of defining basis states and proceeds directly to the de- 
termination of a matrix representation of the algebra of 
observables in which the Hamiltonian is diagonal. The 
approach has its origins in three previous developments: 
a variational technique 0, [?[ for computing the irreps of 
potentially difficult Lie algebras' the double-commutator 
equations-of-motion formalism [J]; and the equations-of- 
motion method of Kerman and Klein 0, Q. Thus, we 
refer to it as the RRKK equations-of-motion method. 

Consider a Hamiltonian H that is a polynomial in the 
elements {X^} of a Lie algebra g of observables with 

commutation relations = J2a^^i^ -^cr- 

fer to as a spectrum generating algebra (SGA), For 
each Xi, e g, define Hi, := [H,X,j]. The objective is to 
determine a unitary irrep in which each observable X^, is 
represented by a matrix X{i'), with elements 



to be determined along with energy differences {Ei~Ej}, 
such that the sets of equations 

f.,{u) := mAX)\j) - {E, - i?,)(z|l.|j) - 0, (3) 

are satisfied. In general, it will also be necessary to in- 
clude additional equations to ensure that the representa- 
tion is the one desired. For example, if the Lie algebra 
has Casimir invariants, C„(g), equations can be included 
to require that they are represented by the appropriate 
multiples of the unit matrix, c„I. 

Let F, a so-called objective function, be defined for a 
finite irrep of dimension N as the sum of squares 

N 



(4) 



X,,{v) := {i\X„\j), 



(1) 



F cannot be negative and can only vanish when a pre- 
cise solution to the system of equations has been ob- 
tained. Thus, for finite irreps, precise solutions to the 
above equations arc obtained by minimization of as a 
function of the unknown matrix elements of the observ- 
ables and the energy differences. If needed, the Hamilto- 
nian matrix can also be evaluated as a polynomial in the 
{X{v)} matrices to determine the ground-state energy. 

The challenge is to obtain accurate solutions for finite 
submatriccs of the observables, corresponding to a sub- 
set of lowest-energy eigenstates, when the irrep is infinite. 
Recall that a differential equation defined over the pos- 
itive half of the real line, < r < oo, can be solved 
precisely over a finite interval < r < i?, if one knows 
the boundary conditions at R. Similarly, when infinite 
dimensional matrices {X(/i)} are truncated to finite di- 
mensions, their outer rows and columns provide bound- 
aries and their entries can be adjusted to give accurate 
results for the matrices they enclose. 

The concern with boundary conditions in the 
equations-of-motion approach arises because, the equa- 
tions of motion for a subset of matrix elements 
{Xij{v); \ < i,j < N}, involve the commutation rcla- 
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tions 



In terms of su(l.l) raising and lowering operators, 



and, hence, Xik{n) matrix elements with k > N. 

This concern can be resolved as follows: because the 
matrix elements connecting lowest-energy states are of 
most interest, we apply a weighting factor Wij = + 
+ 1)] to the expressions that should vanish for an 
exact solution; i.e., redefine the objective function to be 
minimized as 



F' = 



N 



-5]|(»|C„(fl)-c„|j)| 



(6) 



We also make the approximation 

AT+l 

(m, i') « E [^^'^ (^) - ^^fe (m)] (7) 

k=l 

and corresponding approximations for the evaluation of 
fij{v) and the matrix elements of Casimir invariants. 
Minimization of the objective function is then carried 
out iteratively starting from a first guess. Thus, a sim- 
ple physical model of the system can be used to provide 
a (possibly inconsistent) first guess which can then be 
made accurate by the equations-of-motion method. 

To illustrate, consider a Hamiltonian, of relevance in 
the nuclear collective model [§, [l^, [lH , 



Hn 



2M 



f [(l-2a)/32 



(8) 



for an object in a five-dimensional Euclidean space; /3 is 
a radial coordinate in harmonic oscillator units, the 
Laplacian, and M is a dimensionless mass parameter. 
Such a Hamiltonian is of general interest as a model of a 
system with two phases: when a < 0.5, the potential 

Vo.{P) = f[{l-2a)/3' + af3'] (9) 

has a spherical minimum (at /3o = 0); and when a > 0.5, 
it has a minimum given by (Jq ~ {2a — l)/2a. It is 
invariant under the group of SO (5) rotations in the five- 
dimensional space. Thus, its eigenfunctions are products 
of P wave functions and S0(5) spherical harmonics. For a 
state of S0(5) angular momentum v, the /3 wave function 
is an eigenfunction of the radial component of Ha 



Hi-) = -JL_ 



Li! 

2M df3^ 



{v + l)iv + 2) 
2M/32 
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This radial Hamiltonian can be expressed 
of an su(l,l) Lie algebra spanned by 

{v + l){v + 2) 



'2a)f3^+a(3'^]. 

_ (10) 
1121 in terms 



/32 

1 + 2/3 



(11) 



dp 



X 



(v) 



M 



{v) 



MX. 



(v) 



M 



(v) 



(12) 



which satisfy the commutation relations 



■'o ' 



(13) 



we obtain 

Hi"^ = 2(1- 



a 
2M 



(25^") 



S_ 
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[v) 



(14) 



For states of S0(5) angular momentum v, the su(l,l) 
Casimir invariant 



C,{sn{l,l))^So{So-l)-S+S- 
takes the value 



(15) 



(16) 



To start the minimization process, a first guess is 
provided for small a (compared to the critical value 
ttc = 0.5) by the RPA that retains only quadratic terms 
in a Taylor expansion of the Hamiltonian. In the present 
example, this amounts to dropping the quartic term in 
the potential and making the approximation 



Vam^^{l~2a)l3\ 



(17) 



For large a, the Hamiltonian Ha '^ can be approximated 
by its asymptotic limit, obtained from a Taylor expansion 
of VaiP) about its minimum Va(/3o), 



2M dp^ 



{v + l){v + 2) Muj^ 



2MPI 



iP - Po) 



(18) 

This approximation becomes precise as a ^ oo or, for 
a > 0.5, as M oo. The physical content of these two 
limiting solutions is clear. The first is that of a spheri- 
cal harmonic vibrator. The second is that of a rotor in 
a five-dimensional space, with moment of inertia M/3q, 
coupled to a harmonic radial /^-vibrator. Thus, with the 
substitution 



P-Po 



1 



V2Ml 



= (cUc), 



(c-ct). 



d 

dp' \l 2 
the asymptotic Hamiltonian is given by 



Mu) 



(19) 
(20) 

(21) 
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Corresponding approximations for the su(l,l) operators 
are given by 



Ai = ^ (c'c' -f cc — c'c — 



[/3o 



1 



V2M, 



Po 

(ct+c)][;9o+-=i=(ct + c)],(23) 
v2M(' ■ 



X3 = 1 + [\/2M^/3o + (ct + c)] (c - ct). 



(24) 



These approximations, provide first guesses for a < 0.5 — 
d and a > 0.5 + S respectively. In the transition region, 
\a — 0.5| < (S, one can proceed in steps using the results 
of one calculation first guess for the next. 

Before presenting results, it is instructive to consider 
what would be required for accurate results by diagonal- 
ization methods. Fig. [T] shows the expansion coefficients 
for the V = ground state in a harmonic oscillator basis 
for a = 2, 10 and M = 2000. The figure suggests that at 
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FIG. 1: Expansion coefficients for an expansion of the ground 
states of H{a) with a — 2,10 and M = 2000 in the spherical 
(a = 0) harmonic oscillator basis. 

least 420 basis states, for a = 2.0, and 510 basis states, 
for a = 10.0, would be needed for the ground-state wave 
function. More basis states would be required for ex- 
cited states. However, because the equations-of-motion 
approach does not use a predetermined basis, it is pos- 
sible to obtain accurate values for energy eigenvalues for 
low-energy states and matrix elements in a basis of en- 
ergy eigenstates with much smaller matrices. 

The number of unknowns to be determined is re- 
duced by exploiting the fact that S'q"' is self-adjoint and 
(S^^^y = S^^\ The lowest energy level was set to zero 
and the energies {Ei} of excited states regarded as un- 
knowns. For the calculations reported here, the function, 
F' was evaluated using Maple and minimized using the 
'Isqnonlin' algorithm in Matlab. Results obtained with 

= 10 and a range of a and M values are shown be- 
low. Minimization of F' determines 197 unknowns and 
approximately satisfies 500 equations. The values of m 
for which the entries of the m x m submatrices are ac- 
curate to better than 8 significant figures are listed in 
Table HI Numerically precise reference results were ob- 
tained by diagonalizing very large Hamiltonian matrices 
in the spherical harmonic basis. 



TABLE I: Values of m for the m x m submatrices, obtained 
with A'^ = 10 for v = 2, whose elements are accurate to at 
least 8 significant figures. The numbers in brackets are the 
estimated sizes of the Hamiltonian matrix required to achieve 
similar accuracy by diagonalization. 



Q 


Af = 5 


M = 50 


M = 500 


M = 5000 


0.2 


m =7 (25) 


8 (15) 


9 (20) 


9 (20) 


0.5 


7 (30) 


7 (20) 


7(45) 


7(95) 


0.7 


7 (30) 


6 (20) 


6 (75) 


7 (490) 


1.0 


7(35) 


5 (25) 


7 (100) 


8 (770) 


2.0 


7(45) 


5 (35) 


7 (135) 


8 (1065) 


5.0 


6 (65) 


6 (60) 


8 (175) 


8 (1235) 


10.0 


6 (80) 


6 (85) 


8 (200) 


8 (1295) 


50.0 


5 (150) 


7 (180) 


8 (310) 


8 (> 1500) 



Typical computation times to obtain the result shown 
ranged from a few seconds to tens of seconds, achieving 
minimum values of F' ~ 10"^". A satisfying feature of 
the equations-of-motion approach is that its advantages 
over conventional diagonalization are most pronounced 
for large values of M for which the diagonalization ap- 
proach is most slowly convergent. Worst case scenar- 
ios for the equations-of-motion approach are when M is 
small and a is large and when M is large and a = 0.5. 
In the former case, the vibrational fluctuations about 
the equilibrium deformation are large, and in the lat- 
ter case the critical point, a = 0.5, is highly singular. 
Even though the time taken to reach a minimum in such 
situations may be long, the results are invariably accu- 
rate. It is also noteworthy that, in the absense of good 
starting guesses, it is always possible to progress in steps 
from previously found solutions. 

For the present model, it so happens that the a = 0.5 
results are among the easiest to obtain for any value of 
M . This is because of a critical point scaling symmetry 
which means that if the results are known for one 
value of M they can simply be inferred for any M. For 
example, the Hamiltonian at a = 0.5 can be expressed 



1 



i?0.5--^V +-f3 _2_^^^^/3 



1 



+ -/3^ , (25) 



with P = M^^^p. Thus, the energy-level spectrum of 
Ho. 5 is independent of M to within an M~^/^ scale factor. 

Table In] gives an indication of the accuracy obtainable, 
for states of interest, by equations-of-motion calculations. 
It compares excitation energies for the Hamiltonian, Ha 
with a = 2, V ~ 2, and M = 50, obtained for various 
subspace dimensions with those of precise calculations. 
Even a calculation with N = 4 gives accurately the first 
excitation energy to better than one part in 10^. This 
numerical accuracy is sustained for calculated energies 
up to the N — 3 level. 

Table Hill shows the excitation energies obtained by the 
equations-of-motion method with N = 5 for the lowest 
two states of each S0(5) angular momentum v = 0, . . . ,5 
for M = 50 and a range of values of a. The results are 
precise to better than the level of precision shown. 
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TABLE II: Comparison between exact excitation energies and 
values calculated from the equations of motion for a = 2, 
V = 2, and M = 50 for different values of A*'. 



Exact 


N = i 


N = 6 


N = 10 


2.420379968671 
4.799716618277 
7.136170688911 
9.427750466577 
11.672330365421 
13.867713374832 
16.011789014165 
18.102896804066 
20.140619242614 


2.420382 
4.7999 
7.01 


2.420379969 
4.7997167 
7.136180 
9.4280 
11.48 


2.420379968671 
4.799716618277 
7.136170688912 
9.42775046660 
11.672330368 
13.8677136 
16.01180 
18.1033 
20.152 



TABLE III: The lowest two energy levels (n = 0,1) with SO(5) 
angular momentum w = 0, . . . , 5, relative to the ground state 
energy, computed by the equations-of-motion method with 
M = 50 and N — 5. The results agree with those obtained 
by precise calculations to at least the level of precision shown. 
Results of the asymptotic approximation (AS) are shown in 
the last column for a = 10. 



n 


V 


a = 


a = 0.5 


a = 1 


Q = 10 


AS,a = 10 


























1 


1 


0.40223 


0.08712 


0.04254 


0.04211 





2 


2 


0.84184 


0.21629 


0.10634 


0.10526 





3 


3 


1.31343 


0.38597 


0.19140 


0.18947 





4 


4 


1.81319 


0.59436 


0.29769 


0.29474 





5 


5 


2.33829 


0.83960 


0.42521 


0.42105 


1 





2 


0.96133 


1.35973 


6.13301 


6.17262 


1 


1 


3 


1.45802 


1.46666 


6.17646 


6.21472 


1 


2 


4 


1.97813 


1.62318 


6.24161 


6.27788 


1 


3 


5 


2.52043 


1.82581 


6.32847 


6.36209 


1 


4 


6 


3.08360 


2.07109 


6.43701 


6.46735 


1 


5 


7 


3.66642 


2.35579 


6.56723 


6.59367 



It is instructive to note that the domains in which the 
diagonalization and equations-of-motion arc most suc- 
cessful tend to be complementary. Matrix diagonahza- 
tion is often faster. However, it requires much larger 
matrices and substantial extra work to set up the initial 
matrices and interpret the results in a given basis. In 
contrast, the equations-of-motion method, directly com- 
putes the matrices of all observables in the physically 
relevant basis of energy eigenstates. The advantages of 
the equations-of-motion are most evident for states of 
a ^ 0.5 and large M. This is when states are beginning 
to approach their asymptotic dynamical symmetry lim- 
its. For example, for i> = 2, a = 50, and M = 5000, one 
can obtain all the spectroscopic properties of the m = 8 
lowest states in 93 s, using N = 10. Similar accuracy was 
only achieved by diagonalizing larger than 1500 x 1500 
matrices and took 250 s on the same computer. If one 
adds in the time taken to compute the matrices of Hamil- 
tonian and other observables in a predefined basis, the 
equations-of-motion method wins hands down in such a 
situation. These advantages of the RRKK approach are 
expected to be more pronounced for systems with many 
degrees of freedom but relatively simple SGA's. 

An attractive property of the RRKK approach is that 
it makes it possible to start from a simple approxima- 
tion and make it precise. Thus, the RRKK equations 
are particularly relevant for the description of systems 
that exhibit quantum phase transitions with variation of 
parameters, as this letter demonstrates. A natural possi- 
bility is to start with a mean field description of a phase 
transition and, on either side of the critical point, to use 
the RRKK equations of motion to add the fluctuation 
contributions that are missing in the mean field treat- 
ment. 
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